home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / libs / gle / util / surf / hide.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-11-29  |  28.3 KB  |  1,152 lines
  1. /*
  2.     3d Surface plotting with hidden line removal
  3. */
  4. #include "general.h"
  5. #include "gsurface.h"
  6. #include "vdevice.h"
  7. #ifdef unix
  8. #define huge
  9. #endif
  10. #ifdef VAXC
  11. #define huge
  12. #endif
  13. int initminmax(void);
  14. int init_user(void);
  15. int nice_ticks(float *dticks, float *gmin,float *gmax , float *t1,float *tn);
  16. int grid_back(int nx, int ny, float z1, float z2);
  17. int skirt(float huge *z,int ix1, int iy1, float minz);
  18. int find_splits(int nx, int ny, int *splitx, int *splity);
  19. int fxy_polar(float dx,float dy,float *radius,float *angle);
  20. int fpolar_xy(float r, float angle, float *dx, float *dy);
  21. int draw_zaxis(struct axis_struct *ax, int nx, int ny, float minz, float maxz);
  22. int draw_axis(struct axis_struct *ax, int nx, int ny, float minz, float maxz);
  23. int hide(float huge *z,int nx, int ny, float minz, float maxz, struct surface_struct *sff);
  24. int touser3(float x, float y, float z, float *uux, float *uuy, float *uuz);
  25. int touser(float x, float y, float z, float *uux, float *uuy);
  26. int matmove(float i[4][4],float x, float y, float z);
  27. int seth2(int rx1, int ry1, float rz1, int rx2, int ry2, float rz2);
  28. int matscale(float i[4][4],float x, float y, float z);
  29. int matun(float m[4][4]);
  30. int matmul(float a[4][4],float b[4][4]);
  31. int matrx(float i[4][4], float angle);
  32. int matry(float i[4][4], float angle);
  33. int matrz(float i[4][4], float angle);
  34. int vector(int x1, float y1, int x2, float y2);
  35. int draw_riselines(int nx, int ny,float minz, float maxz);
  36. int move3d(float x, float y, float z);
  37. int line3d(float x, float y, float z);
  38. int cube(float x, float y, float z1, float z2);
  39. int horizon(float huge *z,int ix1,int iy1,int ix2,int iy2);
  40. int horizon2(float huge *z,int ix1,int iy1,int ix2,int iy2);
  41. int horizonv(float huge *z,int ix1,int iy1,int ix2,int iy2);
  42. int horizonv2(float huge *z,int ix1,int iy1,int ix2,int iy2);
  43. int hclipvec(int x1, float y1, int x2, float y2, int sethi);
  44. int hclipvec2(int x1, float y1, int x2, float y2, int sethi);
  45. int clipline(float x1, float y1, float z1, float x2, float y2, float z2);
  46. int draw_markers(int nx, int ny);
  47. int matshow(char *s, float m[4][4]);
  48. int show_horizon();
  49. float base;
  50. extern int ngerror;
  51. extern int batch_only;
  52. FILE *fp;
  53. float image[4][4];
  54. init_user()
  55. {
  56.     int i,j;
  57.     matun(image);
  58. }
  59. float smin_x,smax_x,smin_y,smax_y,smin_z,smax_z;
  60. initminmax()
  61. {
  62.     smin_x = 10e10; smax_x = -10e10;
  63.     smin_y = 10e10; smax_y = -10e10;
  64.     smin_z = 10e10; smax_z = -10e10;
  65. }
  66. int setaminmax(float x, float *x1, float *x2);
  67. setaminmax(float x, float *x1, float *x2) 
  68. {
  69.     if (x < *x1) *x1 = x;
  70.     if (x > *x2) *x2 = x;
  71. }
  72. int setminmax(float x,float y,float z);
  73. setminmax(float x,float y,float z)
  74. {
  75.     setaminmax(x,&smin_x,&smax_x);
  76.     setaminmax(y,&smin_y,&smax_y);
  77.     setaminmax(z,&smin_z,&smax_z);
  78. }
  79. #define RERR (0.0001)     /* Rounding error */
  80. int MAXH;
  81. #define deg(d) ((d)*3.14159254/180)
  82. float *h,*h2;     /* h2 is the underneath */
  83. float eye_x,eye_y,eye_depth;
  84. float vdist=0; /* 1=eye touching front of picture, infinity=  cube back */
  85. float map_mul,map_sub;
  86. float maxdepth;
  87. float xmargin = 2;
  88. float ymargin = 1.5;
  89. int doy=true; 
  90. int dox=true;  /* false */
  91. int nnx;
  92. int vsign=1;
  93. #define returng {v_close(); return;}
  94. static struct surface_struct sf;
  95. hide(float huge *z,int nx, int ny, float minz, float maxz, struct surface_struct *sff)
  96. {
  97.     int i,j,x,y,d;
  98.     float r,angle;
  99.     float width,height;
  100.     float ux,uy,y1,y2,step,v,uz,scalex,scaley,ux3,uy3,ux2,uy2;
  101.     static char dummy[90];
  102.     int lasta, isleft, splitx, splity;
  103.     int x1,x2;    
  104.     sf = *sff; /* Make a local copy of all the paramters */
  105.     init_user();
  106.     eye_x = sf.eye_x;
  107.     eye_y = sf.eye_y;
  108.     vdist = sf.vdist;
  109.     base = sf.screenx;
  110.     xmargin = sf.screenx*2/10.0;
  111.     ymargin = sf.screeny*1.5/10.0;
  112.     width = sf.screenx-.5;
  113.     height = sf.screeny-.5;
  114.     if (sf.title!=NULL) height = height -.7;
  115.     dox = sf.xlines_on;
  116.     doy = sf.ylines_on;
  117.     MAXH = sf.maxh;
  118.  
  119.     if (MAXH==0) MAXH = 1000;
  120.     h = malloc(MAXH*sizeof(float));
  121.     h2 = malloc(MAXH*sizeof(float));
  122.     if (h==NULL || h2==NULL) {
  123.         gprint("Was not able to allocate horizon arrays %d \n",MAXH);
  124.         return;
  125.     }
  126.     if (sf.zaxis.min != sf.zaxis.max) {
  127.         minz = sf.zaxis.min;
  128.         maxz = sf.zaxis.max;
  129.     }
  130.     if (ngerror>0) {
  131.         gprint("Press return to continue\n");
  132.         if (batch_only)  gets(dummy); else text_inkey();
  133.         ngerror = 0;
  134.     }
  135.     v_open(sf.screenx,sf.screeny);
  136.     nnx = nx;    /* save dimensions for use in subroutine */
  137.     vsign = 1; /* doing top half */
  138.     for (i=0; i<MAXH ;i++) h[i] = 0.0;  /* zero horizon */
  139.     maxdepth = 0;
  140.     
  141.     /* make it a 10x10x10 cube */
  142.     matmove(image,0.0,0.0,-minz); /* so 0,0,minz = 0,0,0 */
  143.     if (sf.sizez==0) sf.sizez = (sf.sizex+sf.sizey)/2.0;
  144.     matscale(image,sf.sizex/(float) nx,sf.sizey/(float) ny,sf.sizez/(maxz-minz));
  145.             /* positive z comes towards the viewer */
  146.  
  147.     /* 60 50 20 , 80 0 0 */
  148.     matrx(image,deg(sf.xrotate)); /*  clockwise holding right hand side */
  149.     matry(image,deg(sf.yrotate)); /* clockwise holding top */
  150.     matrz(image,deg(sf.zrotate)); /* clockwise holding front */
  151.  
  152.     /* now rotate cube so line from 0,0,0 --> 0,0,1 is virtical on screen */
  153.     touser(0.0,0.0,0.0,&ux,&uy);
  154.     touser(0.0,0.0,1.0,&ux3,&uy3);
  155.     fxy_polar(ux3-ux,uy3-uy,&r,&angle);
  156.     matrz(image,deg(-90.0+angle)); /* clockwise holding front */
  157.         
  158.     /* Normalize image, shift to touch x=0, y=0, z=0 */
  159.     initminmax();
  160.     for (x=0; x<nx; x+=nx-1) {
  161.      for (y=0; y<ny; y+=ny-1) {
  162.         touser3(x,y,minz,&ux,&uy,&uz);
  163.         setminmax(ux,uy,uz);
  164.         touser3(x,y,maxz,&ux,&uy,&uz);
  165.         setminmax(ux,uy,uz);
  166.      }
  167.     }
  168.  
  169.     scalex = 1;
  170.     scaley = 1;
  171.     if (smax_x>width) scalex = width/(smax_x-smin_x);
  172.     if (smax_y>height) scaley = height/(smax_y-smin_y);
  173.     if (scalex<scaley) scaley = scalex;
  174.  
  175.     image[0][3] = image[0][3] - smin_x + xmargin/scaley;
  176.     image[1][3] = image[1][3] - smin_y + ymargin/scaley;
  177.     image[2][3] = image[2][3] - smax_z;
  178.     
  179.         initminmax();
  180.     for (x=0; x<nx; x+=nx-1) {
  181.      for (y=0; y<ny; y+=ny-1) {
  182.         touser3(x,y,minz,&ux,&uy,&uz);
  183.         setminmax(ux,uy,uz);
  184.         touser3(x,y,maxz,&ux,&uy,&uz);
  185.         setminmax(ux,uy,uz);
  186.      }
  187.     }
  188.     
  189.     scalex = 1;
  190.     scaley = 1;
  191.     if (smax_x>width) scalex = width/smax_x;
  192.     if (smax_y>height) scaley = height/smax_y;
  193.     if (scalex<scaley) scaley = scalex;
  194.     matscale(image,scaley,scaley,scaley);    
  195.  
  196.     initminmax();
  197.     for (x=0; x<nx; x+=nx-1) {
  198.      for (y=0; y<ny; y+=ny-1) {
  199.         touser3(x,y,minz,&ux,&uy,&uz);
  200.         setminmax(ux,uy,uz);
  201.         touser3(x,y,maxz,&ux,&uy,&uz);
  202.         setminmax(ux,uy,uz);
  203.      }
  204.     }
  205.  
  206.     
  207.     /* Deepest point will be at   z==-(smax_z-smin_z) */
  208.     maxdepth = smin_z;
  209.     
  210.     v_move(eye_x,eye_y);
  211.     v_set_hei(.3);
  212.  
  213.     /* decide on mapping of ux onto h[] -------------------------- */
  214. #define maph(ux) (((ux)-map_sub)*map_mul)
  215. #define unmaph(ux) (((ux)/map_mul)+map_sub)
  216.     
  217.     x1 = 10000;
  218.     x2 = -10000;
  219.     for (x=0; x<nx; x+=nx-1) {
  220.      for (y=0; y<ny; y+=ny-1) {
  221.          touser(x,y,minz,&ux,&uy);
  222.         if (ux<x1) x1 = ux;
  223.         if (ux>x2) x2 = ux;
  224.          touser(x,y,maxz,&ux,&uy);
  225.         if (ux<x1) x1 = ux;
  226.         if (ux>x2) x2 = ux;
  227.      }
  228.     }
  229.     /* map x1 -- x2  onto 2 -- 998 */
  230.     x1--; x2++;
  231.     map_sub = x1;
  232.     map_mul = (MAXH-100)/((double) x2-x1);
  233.     /* ----------------------------------------------------  */
  234.  
  235.     /* if the closest line of constant X is not one of the edges
  236.        then split object at EYE(X) 
  237.     */
  238.  
  239.     find_splits(nx,ny,&splitx,&splity);
  240.  
  241.     /* if ux of nx,0,0 is less than ux of 0,ny,0 then */
  242.     /* eye is closer to XAXIS then*/ 
  243.     /*   (lines of constant x will be horizontal */
  244.     /*   process lines of constant X in usual order, front to back*/
  245.     /*   After each x line, draw lines of constant y between lastx and nextx */
  246.     /* else, do y first and x betweens.*/
  247.     
  248.     touser(nx,0,0,&ux,&uy); 
  249.     touser(0,ny,0,&ux2,&uy); 
  250.  
  251.  
  252.     /* Now draw bottom half */
  253.     v_color(sf.bot_color);
  254.     v_lstyle(sf.bot_lstyle);
  255.     vsign = -1;    /* reverse tests for bottom half */
  256.     for (i=0; i<MAXH ;i++) h2[i] = 9999.0;  /* zero horizon */
  257.  
  258.     if (sf.bot_on && !sf.skirt_on) {
  259.     touser(nx,0,0,&ux,&uy);
  260.     touser(0,ny,0,&ux2,&uy);
  261.     if (ux < ux2) { 
  262.      for (x=nx-1; x>=0; x--) {
  263.       if (abortkey()) returng;
  264.       if (dox) for (y=splity; y>0; y--) horizonv2(z,x,y,x,y-1); /* x */
  265.       if (dox) for (y=splity+1; y<ny; y++) if (y>0) horizonv2(z,x,y-1,x,y);
  266.       for (y=splity; y>=0; y--) if (x>0) if (doy) horizonv2(z,x,y,x-1,y);
  267.       for (y=splity+1; y<ny; y++) if (x>0) if (doy) horizonv2(z,x,y,x-1,y);
  268.      }        
  269.     } else {
  270.      for (y=0; y<ny; y++) {
  271.       if (abortkey()) returng;
  272.       if (doy) for (x=splitx; x>0; x--) horizonv2(z,x,y,x-1,y); /* y */
  273.       if (doy) for (x=splitx+1; x<nx; x++) if (x>0) horizonv2(z,x-1,y,x,y);
  274.       for (x=splitx; x>=0; x--) if (y<ny-1) if (dox) horizonv2(z,x,y,x,y+1);
  275.       for (x=splitx+1; x<nx; x++) if (y<ny-1) if (dox) horizonv2(z,x,y,x,y+1);
  276.      }        
  277.     }
  278.     } /* if bot_on a&& !skirt_on */
  279.  
  280.  
  281. /* should clip object nicely at zmin,zmax,  do this inside vector */
  282.     vsign = 1;    /* normal tests for top half */
  283.     v_color(sf.top_color);
  284.     v_lstyle(sf.top_lstyle);
  285.     if (sf.top_on) {
  286.     if (ux < ux2) {
  287.      /* Go from x = splitx to zero, then splitx+1 to nx-1 */
  288.      for (x=nx-1; x>=0; x--) {
  289.       if (abortkey()) returng;
  290.       if (dox) for (y=splity; y>0; y--) horizonv(z,x,y,x,y-1); /* X */
  291.       if (dox) for (y=splity+1; y<ny; y++) if (y>0) horizonv(z,x,y-1,x,y); /* Constant X  (this does join, should be earlier )(*/
  292.       for (y=splity; y>=0; y--) if (x>0) if (doy) horizonv(z,x,y,x-1,y); /* Constant Y */
  293.       for (y=splity+1; y<ny; y++) if (x>0) if (doy) horizonv(z,x,y,x-1,y); /* Constant Y */
  294.      }        
  295.     } else {
  296.      /* Go from y = splity to zero, then splity+1 to ny-1 */
  297.      for (y=0; y<ny; y++) { 
  298.       if (abortkey()) returng;
  299.       if (doy) for (x=splitx; x>0; x--) horizonv(z,x,y,x-1,y); /* Constant y */
  300.       if (doy) for (x=splitx+1; x<nx; x++) if (x>0) horizonv(z,x-1,y,x,y); /* Constant y  (this does join, should be first )(*/
  301.       for (x=splitx; x>=0; x--) if (y<ny-1) if (dox) horizonv(z,x,y,x,y+1); /* Constant x */
  302.       for (x=splitx+1; x<nx; x++) if (y<ny-1) if (dox) horizonv(z,x,y,x,y+1); /* Constant x */
  303.      }        
  304.     }
  305.     } /* sf.top_on */
  306.  
  307.  
  308.     v_color(sf.top_color);
  309.     v_lstyle(sf.top_lstyle);
  310.     if (sf.skirt_on) { /* set h2 to bottom of top surface */
  311.         for (x=nx-1; x>0; x--) seth2(x,0,z[x],x-1,0,z[x-1]);
  312.         for (y=0; y<ny-1; y++) seth2(nx-1,y,z[nx-1+y* (long) nx],nx-1,y+1,z[nx-1+(y+1)* (long) nx]);
  313.     }
  314. /*    show_horizon(); */
  315.     if (sf.skirt_on) {
  316.       for (y=splity; y>=0; y--) skirt(z,nx-1,y,minz);
  317.       for (y=splity+1; y<ny; y++) skirt(z,nx-1,y,minz);
  318.       for (x=splitx; x>=0; x--) skirt(z,x,0,minz);
  319.       for (x=splitx+1; x<nx; x++) skirt(z,x,0,minz);
  320.     }
  321.       if (abortkey()) returng;
  322.     if (sf.skirt_on) { /* set h2 to bottom edge of box */
  323.         for (x=nx-1; x>=0; x--) seth2(x,0,minz,x-1,0,minz);
  324.         for (y= -1; y<ny; y++) seth2(nx-1,y,minz,nx-1,y+1,minz);
  325.     }
  326.     /* now define h2[] if bot wasn't drawn */
  327.     if (!sf.bot_on && !sf.skirt_on) {
  328.         for (x=nx-1; x>0; x--) seth2(x,0,z[x],x-1,0,z[x-1]);
  329.         for (y=0; y<ny-1; y++) seth2(nx-1,y,z[nx-1+y* (long) nx],nx-1,y+1,z[nx-1+(y+1)* (long) nx]);
  330.     }
  331.  
  332.     /* Draw the unit cube for reference   -------------------------- */
  333.     if (sf.cube_on) {
  334.         cube(nx-1,ny-1,minz,maxz);
  335.     }
  336.  
  337.     draw_maintitle();
  338.  
  339. /* now lets try and draw the axes */
  340.     draw_axis(&sf.xaxis,nx,ny,minz,maxz);
  341.     draw_axis(&sf.yaxis,nx,ny,minz,maxz);
  342.     draw_zaxis(&sf.zaxis,nx,ny,minz,maxz);
  343.  
  344. /* The back, right and bottom grid lines */
  345.     grid_back(nx,ny,minz,maxz);
  346.  
  347.     draw_markers(nx,ny);
  348.     draw_riselines(nx,ny,minz,maxz);
  349.  
  350.     free(h); free(h2);
  351.     v_close();
  352. }
  353. #define DVAL(a,b) if ((a)==0) a = (b)
  354. draw_maintitle()
  355. {
  356.     v_set_just("BC");
  357.     if (sf.title==NULL) return;
  358.     v_color(sf.title_color);
  359.     DVAL(sf.title_hei,base/30.0);
  360.     v_set_hei(sf.title_hei);
  361.     v_move(sf.screenx/2.0,sf.screeny-sf.title_hei+sf.title_dist);
  362.     v_text(sf.title);
  363. }
  364. draw_axis(struct axis_struct *ax, int nx, int ny, float minz, float maxz)
  365. {
  366.     float ux,uy,ux2,uy2,ux3,uy3,r,a,ta,r2,x,xx,t1,tn;
  367.     char buff[80];
  368.     if (ax->type>1) return;
  369.     if (!ax->on) return;
  370.     if (ax->type==0) {
  371.         touser(0,0,minz,&ux,&uy);
  372.         touser(nx-1,0,minz,&ux2,&uy2);
  373.     } else {
  374.         touser(nx-1,0,minz,&ux,&uy);
  375.         touser(nx-1,ny-1,minz,&ux2,&uy2);
  376.     }
  377.     v_color(ax->color);
  378.     if (!sf.cube_on) {v_move(ux,uy); v_line(ux2,uy2);}
  379.     fxy_polar(ux2-ux,uy2-uy,&r,&a);
  380.     ta = a ;
  381.     a = a - 90;
  382.     if (ax->ticklen == 0) ax->ticklen = base*.001;
  383.     r = ax->ticklen;
  384.     r2 = ax->ticklen+base*.02+ax->dist;
  385.     fpolar_xy(r,a,&ux2,&uy2);
  386.     fpolar_xy(r2,a,&ux3,&uy3);
  387.     
  388.     DVAL(ax->hei,base/60.0);
  389.     v_set_hei(ax->hei);
  390.     v_set_just("TC");
  391.  
  392.     nice_ticks(&ax->step, &ax->min, &ax->max, &t1, &tn);
  393.  
  394.     for (x=t1; x<=.00001+ax->max; x+=ax->step) {
  395.         if (ax->type==0) {
  396.             xx =  (nx-1)*(x-ax->min)/(ax->max - ax->min);
  397.             touser(xx,0,minz,&ux,&uy);
  398.         } else {
  399.             xx =  (ny-1)*(x-ax->min)/(ax->max - ax->min);
  400.             touser(nx-1,xx,minz,&ux,&uy);
  401.         }
  402.         v_move(ux,uy);
  403.         v_line(ux+ux2,uy+uy2);
  404.         v_move(ux+ux3,uy+uy3);
  405.         if (fabs(x)<(.0001*ax->step)) x = 0;
  406.         sprintf(buff,"%g",x);
  407.         v_rotate(ta);
  408.         if (ax->nolast && x> (ax->max - .5*(ax->step))) /*  */ ;
  409.         else if (ax->nofirst && x==t1) /* do nothing */;
  410.         else v_text(buff);
  411.         v_rotate(-ta);
  412.     }
  413.     v_set_just("TC");
  414.     /* Now draw the title for this axis */
  415.     if (ax->title==NULL) return;
  416.     v_color(ax->title_color);
  417.     DVAL(ax->title_hei,base/40.0);
  418.     v_set_hei(ax->title_hei);
  419.     if (ax->type==0) {
  420.         touser((nx-1)/2.0,0,minz,&ux,&uy);
  421.     } else {
  422.         touser(nx-1,(ny-1)/2.0,minz,&ux,&uy);
  423.     }
  424.     DVAL(ax->title_dist,base/17.0);
  425.     r = ax->title_dist;
  426.     fpolar_xy(r,a,&ux2,&uy2);
  427.     v_move(ux+ux2,uy+uy2);
  428.     v_rotate(ta);
  429.     v_text(ax->title);
  430.     v_rotate(-ta);
  431.     
  432. }
  433. draw_zaxis(struct axis_struct *ax, int nx, int ny, float minz, float maxz)
  434. {
  435.     float ux,uy,ux2,uy2,ux3,uy3,r,a,ta,r2,x,xx,t1,tn;
  436.     char buff[80];
  437.  
  438.     if (!ax->on) return;
  439.     touser(0,0,minz,&ux,&uy);
  440.     touser(0,0,maxz,&ux2,&uy2);
  441.     v_color(ax->color);
  442.     if (!sf.cube_on) {v_move(ux,uy); v_line(ux2,uy2);}
  443.     fxy_polar(ux2-ux,uy2-uy,&r,&a);
  444.     ta = a ;
  445.     a = a + 90;
  446.     if (ax->ticklen == 0) ax->ticklen = base*.001;
  447.     r = ax->ticklen;
  448.     r2 = ax->ticklen+base*.02+ax->dist;
  449.     fpolar_xy(r,a,&ux2,&uy2);
  450.     fpolar_xy(r2,a,&ux3,&uy3);
  451.     DVAL(ax->hei,base/60.0);
  452.     v_set_hei(ax->hei);
  453.     v_set_just("RC");
  454.     nice_ticks(&ax->step, &ax->min, &ax->max, &t1, &tn);
  455.     for (x=t1; x<=.0001+ax->max; x+=ax->step) {
  456.         touser(0,0,x,&ux,&uy);
  457.         v_move(ux,uy);
  458.         v_line(ux+ux2,uy+uy2);
  459.         v_move(ux+ux3,uy+uy3);
  460.         if (fabs(x)<(.0001*ax->step)) x = 0;
  461.         sprintf(buff,"%g",x);
  462.         v_text(buff);
  463.     }
  464.     v_set_just("BC");
  465.     /* Now draw the title for this axis */
  466.     if (ax->title==NULL) return;
  467.     v_color(ax->title_color);
  468.     DVAL(ax->title_hei,base/40.0);
  469.     v_set_hei(ax->title_hei);
  470.     touser(0,0,(maxz-minz)/2.0+minz,&ux,&uy);
  471.     DVAL(ax->title_dist,base/17.0);
  472.     r = ax->title_dist;
  473.     fpolar_xy(r,a,&ux2,&uy2);
  474.     v_move(ux+ux2,uy+uy2);
  475.     v_rotate(a-90);
  476.     v_text(ax->title);
  477.     v_rotate(-a+90);
  478. }
  479. move3d(float x, float y, float z)
  480. {
  481.     float ux,uy;
  482.     touser(x,y,z,&ux,&uy);
  483.     v_move(ux,uy);
  484. }
  485. line3d(float x, float y, float z)
  486. {
  487.     float ux,uy;
  488.     touser(x,y,z,&ux,&uy);
  489.     v_line(ux,uy);
  490. }
  491. show_horizon()
  492. {
  493.     int i;
  494.     FILE *f;
  495.     v_color("RED");
  496.     v_move(unmaph(0),h[0]);
  497.     for (i=0;i<900;i++) {
  498.         v_line(unmaph(i),h[i]);
  499.     }
  500.     v_color("BLUE");
  501.     v_move(unmaph(0),h2[0]);
  502.     for (i=0;i<900;i++) {
  503.         v_line(unmaph(i),h2[i]);
  504.     }
  505.  
  506. }
  507. skirt(float huge *z,int ix1, int iy1, float minz)
  508. {
  509.     float ux,uy,v;
  510.     int i,x1,x2;
  511.     float y1,y2,step;
  512.  
  513.     /*
  514.     touser(ix1,iy1,z[ix1+iy1* (long) nnx],&ux,&y1);
  515.     x1 = maph(ux);
  516.     touser(ix1,iy1,minz,&ux,&y2);
  517.     x2 = maph(ux);
  518.     vector(x1,y1,x2,y2);
  519.     */
  520.  
  521.     clipline(ix1,iy1,z[ix1+iy1* (long) nnx],ix1,iy1,minz);
  522. }
  523. horizonv(float huge *z,int ix1,int iy1,int ix2,int iy2)
  524. {
  525.     float ux,uy;
  526.     int x1,x2;
  527.     float y1,y2;
  528.  
  529.     touser(ix1,iy1,z[ix1+iy1* (long) nnx],&ux,&y1);
  530.     x1 = maph(ux);
  531.     touser(ix2,iy2,z[ix2+iy2* (long) nnx],&ux,&y2);
  532.     x2 = maph(ux);
  533.     hclipvec(x1,y1,x2,y2,true);
  534. }
  535. int doclipping=true;
  536. clipline(float x1, float y1, float z1, float x2, float y2, float z2)
  537. {
  538.     float ux,uy,ux2,uy2;
  539.     int ix1,ix2;
  540.  
  541.     touser(x1,y1,z1,&ux,&uy);
  542.     touser(x2,y2,z2,&ux2,&uy2);
  543.     if (!doclipping) { v_move(ux,uy); v_line(ux2,uy2); return;}
  544.     ix1 = maph(ux); ix2 = maph(ux2);
  545.  
  546. /*    printf("hclipvec %d %g  %d %g \n",ix1,uy,ix2,uy2); scr_getch(); */
  547.     if (abs(ix1-ix2)==1) if (fabs(uy2-uy)>.3) ix1 = ix2;
  548.     hclipvec(ix1,uy,ix2,uy2,false);
  549.     hclipvec2(ix1,uy,ix2,uy2,false);
  550. }
  551. horizonv2(float huge *z,int ix1,int iy1,int ix2,int iy2)
  552. {
  553.     float ux,uy;
  554.     int x1,x2;
  555.     float y1,y2;
  556.  
  557.     touser(ix1,iy1,z[ix1+iy1* (long) nnx],&ux,&y1);
  558.     x1 = maph(ux);
  559.     touser(ix2,iy2,z[ix2+iy2* (long) nnx],&ux,&y2);
  560.     x2 = maph(ux);
  561.     hclipvec2(x1,y1,x2,y2,true);
  562. }
  563.  
  564. hclipvec(int x1, float y1, int x2, float y2, int sethi)
  565. {
  566.     float v,sy,step;
  567.     int i,sd,sx,visible;
  568.  
  569.        if (x1==x2) {
  570.     /* printf("VEC %d %g   %d %g h=%g\n",x1,y1,x2,y2,h[x1]); scr_getch();
  571.      */    if (y2<y1) {v = y1; y1 = y2; y2 = v;}
  572.         if (h[x1]<y2) {
  573.             if (h[x1]>y1) y1 = h[x1];
  574.             vector(x1,y1,x2,y2);
  575.             if (sethi) h[x1] = y2;
  576.         }
  577.         return;
  578.        }
  579.        step = (y2-y1)/(x2-x1);
  580.        sd = -1;
  581.        if (x1<x2) sd = 1;
  582.        step = step*sd;
  583.        visible = false;
  584.        for (i=x1, v=y1; sd*i <= sd*x2; i+=sd, v+=step) {
  585.      if (visible) {
  586.         if (h[i]>v) {
  587.         if (sethi) vector(sx,sy,i-sd,v-step);
  588.         else if (fabs(h[i]-v)<.5) { vector(sx,sy,i,h[i]);}    /* i-sd,v) */
  589.         else vector(sx,sy,i-sd,v-step);
  590.         visible = false;
  591.         } else if (sethi) h[i] = v;
  592.      } else {
  593.         if (h[i]<=v+.0001) {
  594.         sx = i; sy = v; visible = true;
  595.         if (!sethi) {if (i!=x1) if (fabs(v-h[i])<.5) sy = h[i];}
  596.         if (sethi) h[i] = v;
  597.         }
  598.      }
  599.        }
  600.        /* draw end part of line */
  601.        if (visible) vector(sx,sy,x2,y2);
  602. }
  603. seth2(int rx1, int ry1, float rz1, int rx2, int ry2, float rz2)
  604. {
  605.     float vx1,vx2,y1,y2;
  606.     int x1,x2;
  607.     float v,sy,step;
  608.     int i,sd,sx,visible;
  609.  
  610.      /*    printf("seth2 %d %d %g, %d %d %g \n",rx1,ry1,rz1,rx2,ry2,rz2);
  611.     scr_getch(); */
  612.     touser(rx1,ry1,rz1,&vx1,&y1);
  613.     touser(rx2,ry2,rz2,&vx2,&y2);
  614.     x1 = maph(vx1); x2 = maph(vx2);
  615.     if (x1<0) x1 = 0;
  616.     if (x2<0) x2 = 0;
  617.     if (x1>MAXH) x1 = MAXH-1;
  618.     if (x2>MAXH) x2 = MAXH-1;
  619.  
  620.        if (x1==x2) {
  621.         if (y2>y1) {v = y1; y1 = y2; y2 = v;}
  622.         if (h2[x1]>y2) h2[x1] = y2;
  623.         return;
  624.        }
  625.        step = (y2-y1)/(x2-x1);
  626.        sd = -1;
  627.        if (x1<x2) sd = 1;
  628.        step = step*sd;
  629.        visible = false;
  630.        for (i=x1, v=y1; sd*i <= sd*x2; i+=sd, v+=step) {
  631.         if (h2[i]>v) h2[i] = v;
  632.        }
  633. }
  634. hclipvec2(int x1, float y1, int x2, float y2, int sethi)
  635. {
  636.     float v,sy,step;
  637.     int i,sd,sx,visible;
  638.  
  639.        if (x1==x2) {
  640. /*    printf("VEC2 %d %g   %d %g h=%g ",x1,y1,x2,y2,h2[x1]); scr_getch();*/
  641.         if (y2>y1) {v = y1; y1 = y2; y2 = v;}
  642.         if (h2[x1]>y2) {
  643.             if (h2[x1]<y1) y1 = h2[x1];
  644.             vector(x1,y1,x2,y2);
  645.             if (sethi) h2[x1] = y2;
  646.         }
  647. /*        scr_getch(); printf("x \n"); scr_getch();                        */
  648.         return;
  649.        }
  650.        step = (y2-y1)/(x2-x1);
  651.        sd = -1;
  652.        if (x1<x2) sd = 1;
  653.        step = step*sd;
  654.        visible = false;
  655.        for (i=x1, v=y1; sd*i <= sd*x2; i+=sd, v+=step) {
  656.         if (visible) {
  657.             if (h2[i]<v) {
  658.                 if (sethi)  vector(sx,sy,i-sd,v-step);
  659.                 else  if (fabs(h2[i]-v)<0.5) { vector(sx,sy,i,h2[i]);}
  660.                 else  vector(sx,sy,i-sd,v-step);
  661.                 visible = false;
  662.             } else if (sethi) h2[i] = v;
  663.         } else {
  664.             if (h2[i]>=v-.0001) {
  665.                 sx = i; sy = v; visible = true;
  666.                 if (!sethi) {if (i!=x1)  if (fabs(v-h2[i])<.5) sy = h2[i];}
  667.                 if (sethi) h2[i] = v;
  668.             }
  669.         }
  670.        }
  671.        /* draw end part of line */
  672.        if (visible) vector(sx,sy,x2,y2);
  673. }
  674. horizon(float huge *z,int ix1,int iy1,int ix2,int iy2)
  675. {
  676.     float ux,uy,v;
  677.     int i,x1,x2;
  678.     float y1,y2,step;
  679.     touser(ix1,iy1,z[ix1+iy1* (long) nnx],&ux,&y1);
  680.     x1 = maph(ux);
  681.     touser(ix2,iy2,z[ix2+iy2* (long) nnx],&ux,&y2);
  682.     x2 = maph(ux);
  683.  
  684.     if (h[x2]-RERR<=y2 && h[x1]-RERR<=y1) {
  685.         vector(x1,y1,x2,y2);
  686.         return;
  687.     }
  688.     if (h[x2]-RERR<=y2 || h[x1]-RERR<=y1) {
  689.       if (h[x1]-RERR>y1) {
  690.         if (x1==x2) {
  691.             vector(x1,h[x1],x2,y2);
  692.             return;
  693.         }
  694.         step = (y2-y1)/(x2-x1);
  695.         if (x1<x2) {
  696.           for (i=x1, v=y1; i<=x2; i++, v+=step) {
  697.             if (h[i]<=v) {
  698.                 vector(i,v,x2,y2);
  699.                 return;
  700.             }
  701.           }
  702.         } else {
  703.           for (i=x1, v=y1; i>=x2; i--, v-=step) {
  704.             if (h[i]<=v) {
  705.                 vector(i,v,x2,y2);
  706.                 return;
  707.             }
  708.           }
  709.         }
  710.       } else {
  711.         if (x1==x2) {
  712.             vector(x1,y1,x2,h[x2]);
  713.             return;
  714.         }
  715.         step = (y2-y1)/(x2-x1);
  716.         if (x1<x2) {
  717.           for (i=x2, v=y2; i>=x1; i--, v-=step) {
  718.             if (h[i]<=v) {
  719.                 vector(x1,y1,i,v);
  720.                 return;
  721.             }
  722.           }
  723.         } else {
  724.           for (i=x2, v=y2; i<=x1; i++, v+=step) {
  725.             if (h[i]<=v) {
  726.                 vector(x1,y1,i,v);
  727.                 return;
  728.             }
  729.           }
  730.         }
  731.       }
  732.     } 
  733. }
  734. horizon2(float huge *z,int ix1,int iy1,int ix2,int iy2)
  735. {
  736.     float ux,uy,v;
  737.     int i,x1,x2;
  738.     float y1,y2,step;
  739.     touser(ix1,iy1,z[ix1+iy1* (long) nnx],&ux,&y1);
  740.     x1 = maph(ux);
  741.     touser(ix2,iy2,z[ix2+iy2* (long) nnx],&ux,&y2);
  742.     x2 = maph(ux);
  743.  
  744.     if (h[x2]+RERR>=y2 && h[x1]+RERR>=y1) {
  745.         vector(x1,y1,x2,y2);
  746.         return;
  747.     }
  748.     if (h[x2]+RERR>=y2 || h[x1]+RERR>=y1) {
  749.       if (h[x1]+RERR<y1) {
  750.         if (x1==x2) {
  751.             vector(x1,h[x1],x2,y2);
  752.             return;
  753.         }
  754.         step = (y2-y1)/(x2-x1);
  755.         if (x1<x2) {
  756.           for (i=x1, v=y1; i<=x2; i++, v+=step) {
  757.             if (h[i]>=v) {
  758.                 vector(i,v,x2,y2);
  759.                 return;
  760.             }
  761.           }
  762.         } else {
  763.           for (i=x1, v=y1; i>=x2; i--, v-=step) {
  764.             if (h[i]>=v) {
  765.                 vector(i,v,x2,y2);
  766.                 return;
  767.             }
  768.           }
  769.         }
  770.       } else {
  771.         if (x1==x2) {
  772.             vector(x1,y1,x2,h[x2]);
  773.             return;
  774.         }
  775.         step = (y2-y1)/(x2-x1);
  776.         if (x1<x2) {
  777.           for (i=x2, v=y2; i>=x1; i--, v-=step) {
  778.             if (h[i]>=v) {
  779.                 vector(x1,y1,i,v);
  780.                 return;
  781.             }
  782.           }
  783.         } else {
  784.           for (i=x2, v=y2; i<=x1; i++, v+=step) {
  785.             if (h[i]>=v) {
  786.                 vector(x1,y1,i,v);
  787.                 return;
  788.             }
  789.           }
  790.         }
  791.       }
  792.     } 
  793. }
  794. vector(int x1, float y1, int x2, float y2)
  795. {
  796.     int i;
  797.     static int savex;
  798.     static float savey;
  799.     float step,v;
  800.  
  801.     if (x2<0 || x1<0) {gprint("Less than zero \n"); exit(1);}
  802.       /*    if (savex!=x1 || savey!=y1) { */
  803.         v_move(unmaph(x1),y1);
  804. /*    }                                       */
  805.     v_line(unmaph(x2),y2);
  806.     savex = x2;
  807.     savey = y2;
  808. }
  809. vectorzz(int x1, float y1, int x2, float y2)
  810. {
  811.     int i;
  812.     static int savex;
  813.     static float savey;
  814.     float step,v;
  815.  
  816.     if (x2<0 || x1<0) {gprint("Less than zero \n"); exit(1);}
  817.     if (x1==x2) {
  818.         if (vsign*h[x1]<vsign*y1) h[x1] = y1;
  819.         if (vsign*h[x2]<vsign*y2) h[x2] = y2;
  820.     } else {
  821.         step = (y2-y1)/(x2-x1);
  822.         if (x1<x2) {
  823.           for (i=x1, v=y1; i<=x2; i++, v+=step) {
  824.             /* printf("! seth i %d v %g  step %g \n",i,v,step); */
  825.             if (vsign*h[i] < vsign*v) h[i] = v;
  826.           }
  827.         } else {
  828.           for (i=x1, v=y1; i>=x2; i--, v-=step) {
  829.             /* printf("! seth i %d v %g  step %g \n",i,v,step); */
  830.             if (vsign*h[i] < vsign*v) h[i] = v;
  831.           }
  832.         }
  833.     }
  834.     if (savex!=x1 || savey!=y1) {
  835.         v_move(unmaph(x1),y1);
  836.     }
  837.     v_line(unmaph(x2),y2);
  838.     savex = x2;
  839.     savey = y2;
  840. }
  841. touser(float x, float y, float z, float *uux, float *uuy)
  842. {
  843.     float uz,ux,uy,p,q;
  844.     ux = image[0][0]*x + image[0][1]*y + image[0][2]*z + image[0][3];
  845.     uy = image[1][0]*x + image[1][1]*y + image[1][2]*z + image[1][3];
  846.     uz = image[2][0]*x + image[2][1]*y + image[2][2]*z + image[2][3];
  847.  
  848.     ux -= eye_x;
  849.     uy -= eye_y;
  850.     if (maxdepth!=0) {
  851.         p = vdist;         /* almost touching infinity */
  852.         q = uz/maxdepth; /* z is negative at deep end , zero at top */
  853.         ux = ux - ux*p*q/(1-p+p*q);
  854.         uy = uy - uy*p*q/(1-p+p*q);
  855. /* (my old transformation)
  856.         ux = uz*(-ux/eye_depth)+ux;
  857.         uy = uz*(-uy/eye_depth)+uy;
  858. */
  859.     }
  860.     *uux = ux + eye_x;
  861.     *uuy = uy + eye_y;
  862. }
  863. touser3(float x, float y, float z, float *uux, float *uuy, float *uuz)
  864. {
  865.     float uz,ux,uy;
  866.     ux = image[0][0]*x + image[0][1]*y + image[0][2]*z + image[0][3];
  867.     uy = image[1][0]*x + image[1][1]*y + image[1][2]*z + image[1][3];
  868.     uz = image[2][0]*x + image[2][1]*y + image[2][2]*z + image[2][3];
  869.  
  870.     ux -= eye_x;
  871.     uy -= eye_y;
  872.     *uux = ux + eye_x;
  873.     *uuy = uy + eye_y;
  874.     *uuz = uz;
  875. }
  876. matscale(float i[4][4],float x, float y, float z)
  877. {
  878.     
  879.     static float c[4][4];
  880.     c[0][0] = x;
  881.     c[1][1] = y;
  882.     c[2][2] = z;
  883.     c[3][3] = 1;
  884.     matmul(i,c);
  885. }
  886. matmove(float i[4][4],float x, float y, float z)
  887. {
  888.     static float c[4][4];
  889.     c[0][0] = 1;
  890.     c[1][1] = 1;
  891.     c[2][2] = 1;
  892.     c[3][3] = 1;
  893.     c[0][3] = x;
  894.     c[1][3] = y;
  895.     c[2][3] = z;
  896.     matmul(i,c);
  897. }
  898. matmul(float a[4][4],float b[4][4])        /* a = a * b */
  899. {
  900.     static float c[4][4],tot;
  901.     int y,xb,x;
  902.     for (y=0;y<4;y++) {
  903.       for (xb=0;xb<4;xb++) {
  904.         tot = 0;
  905.         for (x=0;x<4;x++) tot += a[x][y] * b[xb][x];
  906.         c[xb][y] = tot;
  907.       }
  908.     }
  909.     memcpy(a,c,4*4*sizeof(float));
  910. }
  911. matun(float m[4][4])
  912. {
  913.     int i,j;
  914.     for (i=0;i<4;i++) for (j=0;j<4;j++) m[i][j] = 0;
  915.     for (i=0;i<4;i++) m[i][i] = 1;
  916. }
  917. matrz(float i[4][4], float angle)
  918. {
  919.     /* Rotate around the Z axis */
  920.  
  921.     float m[4][4];
  922.     matun(m);
  923.     m[0][0] = cos(angle);
  924.     m[1][1] = m[0][0];
  925.     m[0][1] = sin(angle);
  926.     m[1][0] = -m[0][1];
  927.     matmul(i,m);
  928. }
  929. matry(float i[4][4], float angle)
  930. {
  931.     /* Rotate around the y axis */
  932.  
  933.     float m[4][4];
  934.     matun(m);
  935.     m[0][0] = cos(angle);
  936.     m[2][2] = m[0][0];
  937.     m[2][0] = sin(angle);
  938.     m[0][2] = -m[2][0];
  939.     matmul(i,m);
  940. }
  941. matrx(float i[4][4], float angle)
  942. {
  943.     /* Rotate around the x axis */
  944.  
  945.     float m[4][4];
  946.     matun(m);
  947.     m[1][1] = cos(angle);
  948.     m[2][2] = m[1][1];
  949.     m[1][2] = sin(angle);
  950.     m[2][1] = -m[1][2];
  951.     matmul(i,m);
  952. }
  953. #define SX(x) ((nx-1)*(x-sf.xaxis.min)/(sf.xaxis.max-sf.xaxis.min))
  954. #define SY(y) ((ny-1)*(y-sf.yaxis.min)/(sf.yaxis.max-sf.yaxis.min))
  955. grid_back(int nx, int ny, float z1, float z2)
  956. {
  957.     float x,y,z;
  958.     v_color(sf.back_color);
  959.     v_lstyle(sf.back_lstyle);
  960.     doclipping = sf.back_hidden;
  961.     if (sf.back_ystep>0) for (y=sf.yaxis.min; y<=sf.yaxis.max+.00001; y+=sf.back_ystep) {
  962.         clipline(0,SY(y),z1,0,SY(y),z2);
  963.     }
  964.     if (sf.back_zstep>0) for (z=z1; z<=z2;  z+=sf.back_zstep) {
  965.       if (abortkey()) return;
  966.         clipline(0,0,z,0,ny-1,z);
  967.     }
  968.  
  969.     v_color(sf.right_color);
  970.     v_lstyle(sf.right_lstyle);
  971.     doclipping = sf.right_hidden;
  972.     if (sf.right_xstep>0) for (x=sf.xaxis.min; x<=sf.xaxis.max+.00001; x+=sf.right_xstep) {
  973.       if (abortkey()) return;
  974.         clipline(SX(x),ny-1,z1,SX(x),ny-1,z2);
  975.     }
  976.     if (sf.right_zstep>0) for (z=z1; z<=z2;  z+=sf.right_zstep) {
  977.       if (abortkey()) return;
  978.         clipline(0,ny-1,z,nx-1,ny-1,z);
  979.     }
  980.  
  981.  
  982.     v_color(sf.base_color);
  983.     v_lstyle(sf.base_lstyle);
  984.     doclipping = sf.base_hidden;
  985.     if (sf.base_xstep>0) for (x=sf.xaxis.min; x<=sf.xaxis.max+.00001; x+=sf.base_xstep) {
  986.         clipline(SX(x),0,z1,SX(x),ny-1,z1);
  987.     }
  988.     if (sf.base_ystep>0) for (y=sf.yaxis.min; y<=sf.yaxis.max+.00001; y+=sf.base_ystep) {
  989.         clipline(0,SY(y),z1,nx-1,SY(y),z1);
  990.     }
  991. }
  992. draw_markers(int nx, int ny)
  993. {
  994.     float *pnt;
  995.     int i;
  996.     pnt = sf.pntxyz;
  997.     if (*sf.marker == 0) return;
  998.     v_color(sf.marker_color);
  999.     DVAL(sf.marker_hei,base/60.0);
  1000.     v_set_hei(sf.marker_hei);
  1001.     for (i=0; i<sf.npnts; i+=3) {
  1002.         move3d(SX(pnt[i]),SY(pnt[i+1]),pnt[i+2]);
  1003.         v_marker(sf.marker);
  1004.     }
  1005. }
  1006. draw_riselines(int nx, int ny,float minz, float maxz)
  1007. {
  1008.     float *pnt=sf.pntxyz;
  1009.     int i;
  1010.  
  1011.     if (sf.riselines!=0) {
  1012.      v_color(sf.riselines_color);
  1013.      v_lstyle(sf.riselines_lstyle);
  1014.      for (i=0; i<sf.npnts; i+=3) {
  1015.         move3d(SX(pnt[i]),SY(pnt[i+1]),pnt[i+2]);
  1016.         line3d(SX(pnt[i]),SY(pnt[i+1]),maxz);
  1017.      }
  1018.     }
  1019.  
  1020.     if (sf.droplines==0) return;
  1021.     v_color(sf.droplines_color);
  1022.     v_lstyle(sf.droplines_lstyle);
  1023.     for (i=0; i<sf.npnts; i+=3) {
  1024.         move3d(SX(pnt[i]),SY(pnt[i+1]),pnt[i+2]);
  1025.         line3d(SX(pnt[i]),SY(pnt[i+1]),minz);
  1026.     }
  1027. }
  1028. cube(float x, float y, float z1, float z2)
  1029. {
  1030.     if (sf.cube_hidden_on) doclipping = true;
  1031.     else doclipping = false;
  1032.     v_color(sf.cube_color);
  1033.     v_lstyle(sf.cube_lstyle);
  1034.  
  1035.     clipline(x,y,z1,0,y,z1);
  1036.     clipline(0,y,z1,0,0,z1);
  1037.  
  1038.     clipline(0,0,z1,0,0,z2);
  1039.     clipline(0,0,z2,0,y,z2);
  1040.     clipline(0,y,z2,0,y,z1);
  1041.     clipline(0,y,z2,x,y,z2);
  1042.     clipline(x,y,z2,x,y,z1);
  1043.  
  1044.     doclipping = false;
  1045.     clipline(0,0,z1,x,0,z1);
  1046.     clipline(x,0,z1,x,y,z1);
  1047.  
  1048.     if (sf.cube_front_on) {
  1049.         clipline(0,0,z2,x,0,z2);
  1050.         clipline(x,0,z2,x,0,z1);
  1051.         clipline(x,0,z2,x,y,z2);
  1052.     }
  1053. /*
  1054.  
  1055.     move3d(0,0,z1);
  1056.     v_color(sf.cube_color);
  1057.     v_lstyle(sf.cube_lstyle);
  1058.     line3d(x,0,z1);
  1059.     line3d(x,y,z1);
  1060.     line3d(0,y,z1);
  1061.     line3d(0,0,z1);
  1062.  
  1063.     line3d(0,0,z2);
  1064.     line3d(0,y,z2);
  1065.     line3d(0,y,z1);
  1066.     move3d(0,y,z2);
  1067.     line3d(x,y,z2);
  1068.     line3d(x,y,z1);
  1069.     if (sf.cube_front_on) {
  1070.         move3d(0,0,z2);
  1071.         line3d(x,0,z2);
  1072.         line3d(x,0,z1);
  1073.         move3d(x,0,z2);
  1074.         line3d(x,y,z2);
  1075.     }
  1076.     */
  1077. }
  1078. matshow(char *s, float m[4][4])
  1079. {
  1080.     int i,j;
  1081.     printf("\n! Matrix {%s} \n",s);
  1082.     for (i=0;i<4;i++)
  1083.         printf("!    %f %f %f %f\n",m[0][i],m[1][i],m[2][i],m[3][i]);
  1084.  
  1085. }
  1086. #define pi (3.141592654)
  1087. find_splits(int nx, int ny, int *splitx, int *splity)
  1088. {
  1089.     int lasta,i,isleft;
  1090.     float ux,uy,ux2,uy2,angle,r;    
  1091.  
  1092.     lasta = 999;
  1093.     *splity = -1;
  1094.     *splitx = nx-1;
  1095.     for (i = 0; i<ny; i++) {
  1096.         touser(nx-1,i,0,&ux,&uy);
  1097.         touser(0,i,0,&ux2,&uy2);
  1098.         fxy_polar(ux2-ux,uy2-uy,&r,&angle);
  1099.         if (angle<90) isleft = true;
  1100.         if (angle>=90) isleft = false;
  1101.         if (lasta==999) lasta = isleft;
  1102.         if (lasta!=isleft) {
  1103.             *splity = i-1;
  1104.         }
  1105.         lasta = isleft;
  1106.     }
  1107.     lasta=999;
  1108.     for (i = 0; i<nx; i++) {
  1109.         touser(i,0,0,&ux,&uy);
  1110.         touser(i,ny-1,0,&ux2,&uy2);
  1111.         fxy_polar(ux2-ux,uy2-uy,&r,&angle);
  1112.         if (angle<90) isleft = true;
  1113.         if (angle>=90) isleft = false;
  1114.         if (lasta==999) lasta = isleft;
  1115.         if (lasta!=isleft) {
  1116.             *splitx = i-1;
  1117.         }
  1118.         lasta = isleft;
  1119.     }
  1120. }
  1121. nice_ticks(float *dticks, float *gmin,float *gmax , float *t1,float *tn)
  1122. {
  1123.     float delta,st,expnt,n;
  1124.     int ni;
  1125.  
  1126.     delta = *gmax-*gmin;
  1127.     if (delta==0) {gprint("Axis range error min=%g max=%g \n",*gmin,*gmax);
  1128.         *gmax = *gmin+10;
  1129.         delta = 10;
  1130.     }
  1131.     st = delta/10;
  1132.     expnt = floor(log10(st));
  1133.     n = st/pow(10,expnt);
  1134.     if (n>5)
  1135.         ni = 10;
  1136.     else if (n>2)
  1137.         ni = 5;
  1138.     else if (n>1)
  1139.         ni = 2;
  1140.     else
  1141.         ni = 1;
  1142.     if (*dticks==0) *dticks = ni * pow(10,expnt);
  1143.     if (*gmin - (delta/1000) <=  floor( *gmin/ *dticks) * *dticks)
  1144.         *t1 = *gmin;
  1145.     else
  1146.         *t1 = (floor(*gmin/ *dticks) * *dticks ) + *dticks;
  1147.  
  1148.     *tn = *gmax;
  1149.     if (( floor( *gmax/ *dticks) * *dticks) < (*gmax - (delta/1000) ))
  1150.         *tn = floor(*gmax/ *dticks ) * *dticks;
  1151. }
  1152.